#pragma once

#include <cstdlib>
#include <cassert>
#include <cmath>
#include <ctime>
#include <initializer_list>
#include <iostream>
#include <utility>
#include <iomanip>

using Real = double;

#define SIGN(a) ((a >= 0) ? 1 : -1)
#define MAX(a, b) ((a) > (b) ? (a) : (b))
#define MIN(a, b) ((a) > (b) ? (b) : (a))

#define ACCURACY 1e-12 // set the ACCURACY
#define ITERATION 100  // set iterate times

/**
 * @file   Matrix.h
 * @author xingyifan <xingyifan@xingyifan-VirtualBox>
 * @date   Wed Jul 14 10:34:03 2021
 *
 * @brief  vector
 *
 *
 */

template <class T>
class Vector
{
protected:
    const int Dim;
    T *value;

public:
    // 初始化构造函数
    Vector(const int Dim, const T &t = T()) : Dim(Dim), value(new T[Dim])
    {
        for (int i = 0; i < Dim; value[i++] = t)
            ;
    }

    // 数组初始化
    Vector(const int Dim, const T arr[]) : Dim(Dim), value(new T[Dim])
    {
        for (int i = 0; i < Dim; value[i] = arr[i], i++)
            ;
    }

    // 可以通过列表初始化
    Vector(std::initializer_list<T> l) : Dim(l.size()), value(new T[Dim])
    {
        auto j = l.begin();
        for (int d = 0; d < Dim; ++d)
            value[d] = *j++;
    }

    // 一定要注意还有普通拷贝构造函数！！！！
    Vector(const Vector &rhs) : Dim(rhs.size()), value(new T[Dim])
    {
        for (int d = 0; d < Dim; d++)
            value[d] = rhs[d];
    }

    // 只能显式地调用此构造函数
    template <class T2>
    explicit Vector(const Vector<T2> &rhs) : Dim(rhs.size()), value(new T[Dim])
    {
        for (int d = 0; d < Dim; d++)
            value[d] = static_cast<T>(rhs[d]);
    }

    // 构造第 D 位为 1 的单位向量
    static Vector<T> unit(int Dim, int D = 0)
    {
        Vector<T> r(Dim);
        r[D] = static_cast<T>(1);
        return r;
    }

    T &operator[](int i) { return value[i]; }
    T operator[](int i) const { return value[i]; }

// 加一个向量
#define VECTOR_BINARY_OP(OpName, Op)             \
    template <class T2>                          \
    auto OpName(const Vector<T2> &rhs) const     \
    {                                            \
        using Tx = decltype(value[0] Op rhs[0]); \
        Vector<Tx> res(Dim);                     \
        for (int i = 0; i < Dim; i++)            \
            res[i] = value[i] Op rhs[i];         \
        return res;                              \
    }

    VECTOR_BINARY_OP(operator+, +)
    VECTOR_BINARY_OP(operator-, -)
    VECTOR_BINARY_OP(operator*, *)
    VECTOR_BINARY_OP(operator/, /)
#undef VECTOR_BINARY_OP

// 加一个常数
#define RIGHT_BROADCAST(OpName, Op)           \
    template <class T2>                       \
    auto OpName(const T2 &rhs) const          \
    {                                         \
        using Tx = decltype(value[0] Op rhs); \
        Vector<Tx> res(Dim);                  \
        for (int d = 0; d < Dim; ++d)         \
            res[d] = value[d] Op rhs;         \
        return res;                           \
    }

    RIGHT_BROADCAST(operator+, +)
    RIGHT_BROADCAST(operator-, -)
    RIGHT_BROADCAST(operator*, *)
    RIGHT_BROADCAST(operator/, /)
#undef RIGHT_BROADCAST

    friend std::ostream &operator<<(std::ostream &out, const Vector<T> &vec)
    {
        for (int i = 0; i < vec.size(); i++)
        {
            out << vec[i] << " ";
        }
        out << std::endl;
        return out;
    }

    void operator=(const Vector &rhs)
    {
        for (int d = 0; d < Dim; d++)
            value[d] = rhs[d];
    }

    // 获取一定范围的向量切片
    Vector<T> range(int left, int right)
    {
        Vector<T> res(right - left);
        for (int i = left; i < right; i++)
        {
            res[i - left] = value[i];
        }
        return res;
    }

    // 检查是否在向量中
    bool find(T t) const
    {
        for (int i = 0; i < Dim; i++)
        {
            if (t == value[i])
            {
                return true;
            }
        }
        return false;
    }

    int size() const { return Dim; }

    ~Vector() { delete[] value; }
};

// ===================================================================================================

// 置换向量（每一位存放所在的位，例如 p[0] = 5，表示 5 元素在 0 位置）
class Perm : public Vector<int>
{
public:
    // 初始化构造函数
    Perm(const int Dim) : Vector<int>(Dim)
    {
        for (int i = -1; i < Dim - 1; this->value[++i] = i)
            ;
    }

    // 拷贝构造
    Perm(const Perm &perm) : Vector<int>(perm.size())
    {
        for (int i = -1; i < Dim - 1; this->value[++i] = perm[i])
            ;
    }

    // 对换
    void exchange(int i, int j)
    {
        int tmp = this->value[i];
        this->value[i] = this->value[j];
        this->value[j] = tmp;
    }

    // 置换的逆
    static Perm inverse(const Perm &p)
    {
        int Dim = p.size();
        Perm q(Dim);
        for (int i = 0; i < Dim; i++)
        {
            q[p[i]] = i;
        }
        return q;
    }

// 置换作用
#define PERM_BINARY_OP(Type, OpName, Op)    \
    auto OpName(const Type &rhs) const      \
    {                                       \
        Type p_rhs(this->Dim);              \
        for (int i = 0; i < this->Dim; i++) \
        {                                   \
            p_rhs[i] = rhs[this->value[i]]; \
        }                                   \
        return p_rhs;                       \
    }

    template <class T>
    PERM_BINARY_OP(Vector<T>, operator*, *)
    PERM_BINARY_OP(Perm, operator*, *)
#undef PERM_BINARY_OP
};

// 如果要在头文件中写函数，最好写成内联函数，避免重复定义
struct Vec
{
    template <class T>
    static bool approx(const Vector<T> &lhs, const Vector<T> &rhs)
    {
        if (rhs.size() != lhs.size())
        {
            return false;
        }
        for (int i = 0; i < lhs.size(); i++)
        {
            if ((lhs[i] - rhs[i]) > ACCURACY)
            {
                return false;
            }
        }
        return true;
    }

    static Vector<Real> rand_vec(int Dim)
    {
        Vector<Real> r(Dim);
        srand((int)time(0));
        for (int i = 0; i < Dim; i++)
        {
            r[i] = rand() / Real(RAND_MAX);
        }
        return r;
    }

    template <class T, int Dim>
    static auto cross(const Vector<T> &va, const Vector<T> &vb)
    {
        if (Dim == 2)
        {
            return va[0] * vb[1] - va[1] * vb[0];
        }
        else if (Dim == 3)
        {
            return Vector<T>{
                va[1] * vb[2] - va[2] * vb[1],
                va[2] * vb[0] - va[0] * vb[2],
                va[0] * vb[1] - va[1] * vb[0]};
        }
        return 0;
    }

    // 幂
    template <class T>
    static Vector<T> power(const Vector<T> &rhs, Real p)
    {
        int Dim = rhs.size();
        Vector<T> res(Dim);
        for (int i = 0; i < Dim; res[i] = pow(static_cast<Real>(rhs[i]), p), i++)
            ;
        return res;
    }

    // 去掉右边的零进行压缩
    template <class T>
    static Vector<T> nonzero(const Vector<T> &rhs)
    {
        int Dim = rhs.size();
        for (int i = Dim - 1; i >= 0; i--)
        {
            if (rhs[i] != 0)
            {
                return rhs.range(0, i + 1);
            }
        }
        return rhs;
    }

    template <class T>
    static Vector<T> reverse(const Vector<T> &rhs)
    {
        int Dim = rhs.size();
        Vector<T> res(Dim);
        for (int i = 0; i < Dim; res[i] = rhs[Dim - i - 1], i++)
            ;
        return res;
    }

    template <class T>
    static T sum(const Vector<T> &rhs)
    {
        T a = 0;
        for (int i = 0; i < rhs.size(); i++)
        {
            a += rhs[i];
        }
        return a;
    }

    template <class T>
    static T dot(const Vector<T> &vec1, const Vector<T> &vec2, Real weight = 1)
    {
        int Dim = vec1.size();
        T sum = 0;
        for (int i = 0; i < Dim; i++)
        {
            sum += vec1[i] * vec2[i] * weight;
        }
        return sum;
    }

    template <class T>
    static Vector<int> sgn(const Vector<T> &rhs)
    {
        int Dim = rhs.size();
        Vector<int> res(Dim);
        for (int i = 0; i < Dim; i++)
        {
            res[i] = (rhs[i] > 0 ? 1 : (rhs[i] < 0 ? -1 : 0));
        }
        return res;
    }

    template <class T>
    static Vector<int> abs(const Vector<T> &rhs)
    {
        int Dim = rhs.size();
        Vector<int> res(Dim);
        for (int d = 0; d < Dim; res[d] = fabs(rhs[d]), ++d)
            ;
        return res;
    }

    template <class T>
    static T norm(const Vector<T> &lhs, int nt = 2)
    {
        int Dim = lhs.size();
        if (nt == 2)
        {
            return sqrt(dot(lhs, lhs));
        }
        else if (nt == 1)
        {
            T a = 0;
            for (int i = 0; i < Dim; a += fabs(lhs[i]), ++i)
                ;
            return a;
        }
        else if (nt == 0)
        {
            T a = fabs(lhs[0]);
            for (int i = 1; i < Dim; i++)
            {
                a = MAX(fabs(lhs[i]), a);
            }
            return a;
        }
        return 0; // never reached
    }

    template <class T>
    static Vector<T> normalize(const Vector<T> &rhs, int nt = 2)
    {
        T normalizer = norm(rhs, nt);
        return Vector<T>(rhs / normalizer);
    }

    template <class T, class Cmp = std::less<T>>
    static T min_of(const Vector<T> &v, Cmp cmp = Cmp())
    {
        int Dim = v.size();
        auto a = v[0];
        for (int d = 1; d < Dim; d++)
            if (cmp(v[d], a))
                a = v[d];
        return a;
    }

    template <class T>
    static T max_of(const Vector<T> &v)
    {
        return min_of(v, std::greater<T>());
    }

    template <class T>
    static int max_index_of(const Vector<T> &v)
    {
        int Dim = v.size();
        auto a = v[0];
        int index = 0;
        for (int d = 1; d < Dim; d++)
        {
            index = (fabs(v[d]) > fabs(v[index]) ? d : index);
        }
        return index;
    }

    // element-wise min
    template <class T>
    static Vector<T> min(const Vector<T> &lhs, const Vector<T> &rhs)
    {
        int Dim = lhs.size();
        Vector<T> res(Dim);
        for (int d = 0; d < Dim; d++)
        {
            res[d] = MIN(lhs[d], rhs[d]);
        }
        return res;
    }

    // element-wise max
    template <class T>
    static Vector<T> max(const Vector<T> &lhs, const Vector<T> &rhs)
    {
        int Dim = lhs.size();
        Vector<T> res(Dim);
        for (int d = 0; d < Dim; d++)
        {
            res[d] = MAX(lhs[d], rhs[d]);
        }
        return res;
    }
};